# Model 1, no covariate adjustment
reg_res_mod1 <- function(outcome="", dataIn=NULL){
  Model1 <- lm(as.formula(paste0(outcome,"~","treat+as.factor(year)")),data=dataIn)
  results <- tidy(Model1)
  results <- results[results[,"term"]=="treat",]
  results$estimate <- coeftest(Model1 , vcov = vcovHC(Model1, type = "HC1"))[2,1]/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$std.error <- coeftest(Model1 , vcov = vcovHC(Model1, type = "HC1"))[2,2]/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$N <- nrow(dataIn)
  results$outcome <- paste0(outcome)
  return(results[c("estimate","std.error","p.value", "N","outcome")])
}  

reg_res_mod1y <- function(outcome="", dataIn=NULL){
  Model1 <- lm(as.formula(paste0(outcome,"~","treat")),data=dataIn)
  results <- tidy(Model1)
  results <- results[results[,"term"]=="treat",]
  results$estimate <- coeftest(Model1 , vcov = vcovHC(Model1, type = "HC1"))[2,1]/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$std.error <- coeftest(Model1 , vcov = vcovHC(Model1, type = "HC1"))[2,2]/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$N <- nrow(dataIn)
  results$outcome <- paste0(outcome)
  return(results[c("estimate","std.error","p.value", "N","outcome")])
}  



# Model 2, controlling for the frequency of weekly FB usage
reg_res_mod2<- function(outcome="", dataIn=NULL){
  covIn <- as.formula("~as.factor(year)+fb_weeklyusage")
  Model2 <- lm_lin(as.formula(paste0(outcome,"~","treat")), covariates=covIn, data=dataIn)
  results <- tidy(Model2)
  results <- results[results[,"term"]=="treat",]
  results$estimate2 <- results$estimate/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$std.error2 <- results$std.error/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$N <- Model2$nobs
  return(results[c("estimate2","std.error2","p.value", "N","outcome")])
}  


# Model 2, controlling for the frequency of weekly FB usage, without the year control
reg_res_mod2y<- function(outcome="", dataIn=NULL){
  covIn <- as.formula("~fb_weeklyusage+age")
  Model2 <- lm_lin(as.formula(paste0(outcome,"~","treat")), covariates=covIn, data=dataIn)
  results <- tidy(Model2)
  results <- results[results[,"term"]=="treat",]
  results$estimate2 <- results$estimate/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$std.error2 <- results$std.error/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$N <- Model2$nobs
  return(results[c("estimate2","std.error2","p.value", "N","outcome")])
}  


# Model 3, controlling for a rich set of covariates
reg_res_mod3 <- function(outcome="", dataIn=NULL){
  covIn <- as.formula("~fb_weeklyusage+gender+age+employ+imp_ethn+imp_cntry+as.factor(year)")
  Model3 <- lm_lin(as.formula(paste0(outcome,"~","treat")), covariates=covIn, data=dataIn)
  results <- tidy(Model3)
  results <- results[results[,"term"]=="treat",]
  results$estimate3<-results$estimate/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$std.error3<-results$std.error/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$N <- Model3$nobs
  return(results[c("estimate3","std.error3","p.value", "N","outcome")])
}  

# Model 3, controlling for a rich set of covariates
reg_res_mod3y <- function(outcome="", dataIn=NULL){
  covIn <- as.formula("~fb_weeklyusage+gender+age+employ+imp_ethn+imp_cntry")
  Model3 <- lm_lin(as.formula(paste0(outcome,"~","treat")), covariates=covIn, data=dataIn)
  results <- tidy(Model3)
  results <- results[results[,"term"]=="treat",]
  results$estimate3<-results$estimate/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$std.error3<-results$std.error/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$N <- Model3$nobs
  return(results[c("estimate3","std.error3","p.value", "N","outcome")])
}  



# Model 3, controlling for a rich set of covariates
reg_res_mod3_noyear <- function(outcome="", dataIn=NULL){
  covIn <- as.formula("~fb_weeklyusage+gender+age+employ+imp_ethn+imp_cntry")
  Model3 <- lm_lin(as.formula(paste0(outcome,"~","treat")), covariates=covIn, data=dataIn)
  results <- tidy(Model3)
  results <- results[results[,"term"]=="treat",]
  results$estimate3<-results$estimate/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$std.error3<-results$std.error/sd(dataIn[dataIn$treat==0,outcome], na.rm=TRUE)
  results$N <- Model3$nobs
  return(results[c("estimate3","std.error3","p.value", "N","outcome")])
}  


# Find the modal value
Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

# PCA Index
#pc_index <- function(dataIn=NULL,
                    # varList=NULL){
 # x <- dataIn[,varList]
 # res <- principal(x, nfactors = 1)
  #index <- res$scores[,1]
  #return(index)
#}

pc_index <- function(dataIn=NULL,
                           varList=NULL){
  x <- dataIn[,varList]
  res   <- princomp(na.omit(x))
  index <- as.vector(predict(res, newdata=x)[,1])
  return(index)
}


# Impute group mode
imputeGroupMode <- function(dataIn = rctd, 
                            varIn = NULL,
                            groupVar = NULL){
  index <- 1:nrow(dataIn)
  varUp <- dataIn[,varIn]
  NA_index <- index[is.na(varUp)]
  for(j in 1:length(NA_index)){
    groupUp <- dataIn[NA_index[j], groupVar]
    varUp[NA_index[j]] <- Mode(na.omit(dataIn[dataIn[,
                                                     groupVar]==groupUp,
                                              varIn]))
  }
  return(varUp)
}

